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ABSTRACT 

The solar and extra solar gas giants appear to have diverse internal structure and metallicities. 
We examine a potential cause for these dispersions in the context of the conventional sequential 
accretion formation scenario, which assumes the initial formation of cores from protoplanetary 
embryos. In principle, gas accretion onto cores with masses below several times that of the Earth 
is suppressed by the energy released from the bombardment of residual planetesimals. After the 
cores have attained their isolation masses, additional mass gain through gas accretion enlarges 
their feeding zones and brings a fresh supply of planetesimals. However, the relatively low-mass 
cores have limited influence on exciting the eccentricities of the newly embraced planetesimals. 
Due to their aerodynamical and tidal interaction with the nascent gas disk, planetesimals on 
eccentric orbits undergo slow orbital decay. We show that these planetesimals generally cannot 
pass through the mean motion resonances of the cores, and the suppression of planetesimal bom- 
bardment rate enables the cores to accrete gas with little interruption. During growth from the 
cores to protoplanets, the resonance width of protoplanets increases with their masses. When 
the resonances overlap with each other, the trapped planetesimals become dynamically unstable 
and their eccentricity excitation is strongly enhanced. Subsequent gas drag induces the plan- 
etesimals to migrate to the proximity of the protoplanets and collide with them. This process 
leads to the resumption and a surge of planetesimal bombardment during the advanced stage of 
the protoplanet growth. The most massive intruders are the residual earth-mass protoplanetary 
embryos that may survive the passage through the protoplanet envelopes and increase their core 
masses. This mechanism may account for the diversity of the core-envelope structure between 
Jupiter, Saturn and the metallicity dispersion inferred from the transiting extra solar planets. 
During the final formation stage of the proto-gas-giants, their tidal torque induces the formation 
of gaps in the gas disk. This perturbed structure of gas disk also leads to the accumulation of 
planetesimals outside the feeding zone of the protoplanets. The surface density enhancement 
promotes the subsequent buildup of cores for secondary gas giant planets outside the orbit of the 
first-born protoplanets and the formation of eccentric multiple planet systems. 

Subject headings: Giant planet; planetary formation; solar nebula; planet dynamics; n-body simulation 



1. Introduction 

Models of the interior structure of giant plan- 
ets in our Solar System suggest that giant planets 
contain heavy elements with total masses up to 
several tens Earth masses (M^) (e.g., Wuchterl 
et al. 2000). Due to uncertainties in the equa- 
tion of state, the internal distribution of the heavy 
elements is poorly determined. Guillot, Gautier 



& Hubbard (1997) constructed models of Jupiter 
which indicate the presence of a core with mass in 
the range — 12M{^ and a total mass of heavy ele- 
ments reaching 11— 45M®. With a slight modifica- 
tion in the equation of state used in these models, 
Saumon & Guillot (2004) deduced a smaller up- 
per limit for the core mass of Jupiter. They find a 
Jupiter model with core mass of 0— IIM^ and to- 
tal heavy elements of 8— 39M0. Their models sug- 
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gest Saturn may have a core with mass 9 — 22M0 
and a total heavy-element mass of 13 — 28Mq. 

Super-solar metallicity and massive core have 
also been inferred from the transit observation 
of a compact Saturn- mass extra solar planet, 
HD149026b (Sato et al. 2005). However, transit 
observations of several other close-in extra solar 
planets indicate a large dispersion in the plane- 
tary mass-size distribution. Since the contraction 
rate during the evolution of gas giant planets is 
determined by the ratios of their cores to envelope 
masses (Bodenheimer et al. 2001, Burrows et al. 
2000), a large spread in their mass-radius relation 
is indicative of a wide dispersion in their internal 
structure. 

The existence of massive cores in Saturn and 

HD149G26b provides a strong support for the con- 
ventional sequential accretion scenario (Perri & 
Cameron 1974, Mizuno 1980, Bodenheimer & Pol- 
lack 1986, Pollack et al. 1996), which is based 
on the assumption that gas giant planets form 
through three major stages: 

51) Embryo-growth stage: protoplanctary cores 
formed and grew mainly by the bombardment 
of planetesimals onto them. Although low-mass 
cores (up to a few M@) can also attract a small 
amount of gas, the envelopes initially built up at 
a much slower rate than the growth of the core. 
The cores attain isolation masses at the end of 
this initial stage. 

52) Quasi-hydrostatic sedimentation stage: the 
accretion of planetesimals tapers as their supply in 
the feeding zone is depleted by the cores. As the 
dissipation rate also decreases along with the de- 
clining influx of planetesimal bombardment, ther- 
mal energy continues to diffuse out of the envelope. 
This loss of entropy allows a quasi- hydrostatic sed- 
imentation and the growth of the gaseous enve- 
lope. 

53) Runaway gas-accretion stage: when the mass 
of the gas becomes comparable to that of the core, 
the rate of gas sedimentation increases with the 
intensified flux of radiation transfer. The char- 
acteristic growth time scale of the protoplanets 
decreases with their masses. This runaway stage 
continues until the gas supply is exhausted by ei- 
ther the formation of a tidally induced gap near 
the protoplanet orbit or the depletion of the entire 
nascent disk. 



While this paradigm has been widely accepted, 
many uncertainties remain. One major issue con- 
cerns the protracted transition through the sec- 
ond stage (S2). In the early models (Pollack et al. 
1996), this stage persists over a time scale longer 
than the observationally inferred depletion time 
scale of the disk (^ a few Myr), because any in- 
crease in the protoplanet mass also leads to an 
expansion of its feeding zone and a resurgence of 
planetesimal accretion, which tends to slow down 
the gas accretion. A simple extrapolation of the 
early models would imply a low probability of gas- 
giant formation, which is incompatible with the 
observationally inferred ubiquity (with a proba- 
bility ~ 0.1 — 0.15) of gas giant planets around 
nearby solar- type stars (Marcy 2005). 

Another unsolved issue is the compositional di- 
versity. Regardless of the largo dispersion in the 
core mass, the average metallicity of Jupiter is 
about twice solar while that of Saturn is an order 
of magnitude larger than that of the Sun. There 
are several potential mechanisms which may lead 
to these differences: 

Ml) The critical mass of the cores (Mcr) needed 
for the onset of efficient gas accretion (i.e. from 
quasi-hydrostatic sedimentation (S2) to runaway 
gas-accretion stage (S3)) depends on their plan- 
etesimal bombardment rate (Mp) and the rate of 
radiation transfer (Ikoma et al. 2000). Both pro- 
cesses may be stochastic and dependent on the 
inventory of residual planetesimals and the dust- 
to-gas ratio. 

M2) During the runaway gas-accretion stage (S3), 
a fraction of the pre-existing cores may be eroded 
and mixed into the envelope. 
M3) The metallicity of the accreted gas in the run- 
away gas-accretion stage (S3) may depend on the 
epoch of protoplanet formation, since its value in- 
creases during the transition from protostellar to 
debris disks. 

M4) During and after the runaway gas-accretion 
stage (S3), gaseous planets may also gain mass in 
heavy elements through giant impacts by nearby 
residual planetesimals and protoplanctary em- 
bryos. Most intruders are disrupted during their 
passage through the envelope. However, colliding 
embryos with several may be able to reach 
the cores and increase their masses. 

Here, we explore how several relevant physical 
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processes may act together to overcome the growth 
challenge of gas giants and introduce metallicity 
diversity. The main aim of this paper is to examine 
the dynamical interaction of a growing proto-gas- 
giant planet with its neighboring planetesimals 
and protoplanetary embryos in a gaseous environ- 
ment. This analysis is relevant in two contexts. 
Our first objective is to evaluate whether process 
Ml can significantly reduce the planetesimals ac- 
cretion rate onto relatively low-mass (a few Mg) 
cores. A reduction of the energy dissipation associ- 
ated with planetesimal accretion would enable the 
gas to settle onto the cores at a more rapid speed 
than the early models, thus shorten the transi- 
tion from stages of quasi-hydrostatic sedimenta- 
tion (S2) to nmaway gas-accretion (S3). The re- 
duced influx of the planetesimals would also lower 
the replenishment of dust and the contribution to 
the opacity in the envelope. Protoplanetary mod- 
els show that the suppression of opacity enhances 
both the heat flux through the radiative region 
and the gas accretion rate (Ikoma, et al. 2000, 
Hubickyj et al. 2005). 

This promising avenue to bypass the gas- 
accretion barrier also implies a limited acquisi- 
tion of heavy elements during the initial evolution 
of proto planets when their mass Mp is only a 
few Earth masses. In order to account for the 
rich abundance of heavy elements in the gas gi- 
ants, especially in their envelopes, we need to 
consider a possible mechanism for the resump- 
tion of planetesimal-accretion during the runaway 
gas-accretion stage (S3) when the protoplanets 
have acquired a major fraction of their present- 
day gaseous envelopes (with Mp ~ Mj where Mj 
is Jupiter's mass). Our second goal is to assess 
the efficiency of planetesimal accretion under the 
influence of process M4. The first objective has 
implications for the ubiquity of gas giants whereas 
the second is linked to the structural diversity of 
gas giants. 

In §2 of this paper, we first present a dynami- 
cal model for this process and estimate the various 
time scales associated with both aerodynamical 
and tidally-induced gas drag. In order to simplify 
the problem, we assume a prescribed model for the 
evolution of protoplanet mass which is based on 
the Bondi formula for idealized, spherically sym- 
metric, unimpeded accretion. We also utilize an 
ad hoc uniform accretion prescription to illustrate 



the dominant physical processes which determine 
the dynamical evolution of the system. In §3, 
we show that the combined effects of the proto- 
planct's perturbation and the planetesimals' aero- 
dynamical and tidal interaction with their nascent 
disks lead to the formation of a planetesimal gap 
when Mp is a few Af®. With a reduced rate of 
planetesimal bombardment at the onset of the gas 
accretion, the envelopes contract more rapidly and 
the gas accretion time scale can be shortened from 
the early models. In §4, we show that the plan- 
etesimal accretion rate increases rapidly at the late 
stage of the protoplanct's gas accretion. We also 
provide evidence to demonstrate the accumula- 
tion of planetesimals outside the asymptotic feed- 
ing zone. This effect can promote the formation 
of second- generation proto-gas-giant planets. Fi- 
nally, we summarize our results and discuss their 
implications in §5. 

2. Dynamical model 

In this section, we study the growth of a pro- 
toplanet from the beginning of quasi-hydrostatic 
sedimentation stage (S2), i.e., after its core reaches 
an isolation mass. The protoplanet is moving in a 
gaseous and planetesimal disk. The orbits of the 
protoplanet and planetesimals are perturbed by 
the gas drag. We first estimate the effects of gas 
drag in §2.1. Models of the protoplanet and plan- 
etesimals are given in §2.2 and §2.3, respectively. 

2.1. Gravitational gas drag 

There arc three physical processes that are act- 
ing on the protoplanet and planetesimals during 
their evolution: aerodynamical gas drag, gravita- 
tional tidal drag and dynamical friction. Dynami- 
cal friction on the most massive embryos by the 
low-mass planetesimals plays an important role 
during the stage of runaway growth of protoplan- 
etary embryos (e.g., Wetherill & Stewart 1989, 
Palmer et al. 1993, Kokubo & Ida 1996, Goldre- 
ich ct al. 2004). During the subsequent oligarchic 
stage, massive embryos emerge to perturb the ve- 
locity dispersion (a) of the residual planetesimals 
(Ida & Makino 1993, Kokubo & Ida 1998). Nu- 
merical simulations show that the collisions gen- 
erally lead to coagulation and the embryos attain 
most of the mass in heavy elements (Kokubo & Ida 
2000, Leinhardt & Richardson 2005). In this limit. 
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we assume that the effect of dynamical friction can 
be incorporated into the gravitational drag due to 
embryo-gas interaction. 

2.1.1. Aerodynamical gas drag 

The gas drag on small particles is in the form 
of aerodynamical drag (e.g., Adachi et al. 1976, 
Tanaka & Ida 1999). The acceleration of a plan- 
etesimal with mass m by the aerodynamical drag 
has the form 



Facro = -T^Cz37r5V<,|U|U, 

2m 



(1) 



where Co — 0.5 is the drag coefficient for objects 
with large Reynold's number, S is the radius of the 
planetesimal, pg is the density of gas, U = Vk— Vg 
is the relative velocity, and Vg are the velocity 
vectors of the planetesimal's Keplerian motion and 
gas motion respectively. 

The motion of the gas is subject to both the 
stellar gravity and its own pressure gradient. In 
an unperturbed (by the protoplanet) region of the 
disk, the balance of forces in the radial direction 
gives 

1 dP 



y2 
R 



^ c 

R 



(2) 



where Vc — \JGM^IR is the pressure-free circular 
velocity at the radial distance R to the host star, 
and P is the gas pressure. In a stable disk, P = 
PgCg, where Cs = (^^m^r)^^^ ^'^^ sound speed of 
an ideal gas in an isothermal environment, with fc 
the Boltzmann constant, mu the proton mass, T 
the gas temperature and /z the average molecular 
weight of the gas. We take /x = 2.34 for the gas 
with solar composition, thus Cs = 5.95 x 10'^ cm 
s"V^/K- Suppose that 



„2 _ „2 ( R \-Sr 
~ ^sO \ lAU 



M, 

Mq 



(3) 



where the subscript o denotes the value of the cor- 
responding quantity at lAU (and = Mq in the 
equation of T). From equation ^ we get 



where 



V^ = V,{l-2^{R)Y'\ 



{Sp + 5c) / Cs ^2 

2 ' 



(4) 



(5) 



and 

Vj \Kj \1A\JJ \Mq 

(6) 

Throughout this paper, we use the minimum 
mass nebula model (Hayashi et al. 1985) to eval- 
uate fiducial model parameters, though our basic 
algorithm can be applied to a more general disk. 
In this model, the gas is heated to an equilibrium 
temperature by the central star with temperature 



-1/2 



(7) 



The surface density of the gas disk is given by 
/ R \ ^'^/^ 



(8) 



where E^o = 1700 g cm~^, fg is a scaling factor 
so that /g = 1 corresponds to the minimum mass 
nebula. The corresponding density of gas disk is 



Pg = 1.4 X lO-^g cm-3/. 



'n lAU, 



(9) 



Since we have Sp — 11/4, Sc = 1/2, substituting 
them into equations ([5]) and ([6]) , rj has the form of 



r]{R) = 0.0018 



R 
lAU 



1/2 



(10) 



In more realistic models, the magnitude of rj can 
be modified by the surface irradiation, internal 
viscous dissipation, and the radiation transfer 
through the disk (Garaud & Lin 2007). 

Assuming the protoplanet or a planetesimal has 
a density p, its radius can be expressed as 



M \ ^/ V 1 K cm"3 



1/3 



S* = 5.2 X lO^^AU — 

(11) 

In terms of the above expression, equation ([T]) can 
be expressed as 

/ .r \ -1/3 / \ -2/3 

{^) (rr^T-) Iu|u. 

(12) 

The aerodynamical gas drag decreases the semi- 
major axis a, eccentricity e, and inclination i of the 
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planetesimal orbit. The average time scales for 
the evolution of these orbital elements are given 
by Adachi et al. (1976) 



1 _ 1 tda^ 
^a,a a \dt / aero 



1 _ 1 ( de\ _ 2 r di\ 

e \ dt / aero i \dt / aero 



where Tacro is a time scale given as: 



(13) 



7rC_DS^/3gVk(a) 



3.5x10 



(14) 

Note that Ta.a ^ Ta,e/{V + 6^ + i^) >> Ta,e- 

2.1.2. Tidal effect of gas disk 

Tidal interaction between a gas disk and a pro- 
toplanet leads to an effect similar to gas drag, 
which is particularly important for the dynamical 
evolution of protoplanets with large masses (Gol- 
dreich & Tremaine 1980, Ward 1986, 1989, 1997, 
Artymowicz 1993). We adopt the form of acceler- 
ation from Kominami and Ida (2002): 



2/3 



13/4 



ftidal — 



Vfe - v„ 

Ttidal 



(15) 



where Vg is the velocity of gas motion, and we 
consider it in circular orbits, and rtidai is the time 
scale defined as (Ward 1989, Artymowicz 1993) 



Ttidal 



SxlO'yr 



M 



(ft); 



(16) 



where VL^ is the Keplerian frequency of circular 
motion. Since dynamic friction of planetesimals 
on protoplanets have similar expressions of accel- 
eration and time scale as those of the disk tide, i.e., 
equations p5|) and (fTBj) (see Appendix of Komi- 
nami & Ida 2002), we do not consider the effect of 
dynamical friction particularly in this work. 

Tidal drag decreases the eccentricities and in- 
clinations of the embryo orbits. In principle, tidal 
interaction between embryos and the disk also 
lead directly to the decay of embryo orbits (Ward 
1997). The rate of this "type I migration" is de- 
termined by an imbalance in the torque from disk 



regions interior and exterior to the embryo orbits. 
Linear analysis for this process is evaluated for 
idealized background surface density (Tanaka et 
al. 2002). The results of linear analysis imply 
that, in a solar nebula environment, embryos more 
massive than the Earth would migrate rapidly to- 
ward their host stars, thus gas giants would rarely 
form (Ida & Lin 2007). In general, both intrinsic 
(Rice & Armitage 2003, Laughlin et al. 2004, Nel- 
son & Papaloizu 2003) and self-excited turbulence 
(KoUer et al. 2003) may lead to nonlinear evolu- 
tion of disk structure, readjustment of the torque 
distribution, and the suppression of type I migra- 
tion. However, these structural adjustments do 
not modify the efficiency of eccentricity damping 
since the contribution from both sides of the disk 
is cumulative rather than cancelling. Thus, we in- 
clude here the effect of tidally-induced eccentric- 
ity damping but neglect that of type I migration. 
Nevertheless, there is an associated change in the 
semi- major axes of the planetesimals and embryos. 
Within 0{e^,i^), the average time scales for the 
evolution of these orbital elements are given by. 



1 — 1 fda\ 
'''t.a a \ dt ) tidal 

1 — 1 rde\ 
Tt.e e \dt / tidal 



S'^'tidal 
1 / 



(5e2 + 2i^) 



1 



13„2 
32 



1 ldA\ 
i Vdt/ti 



dal 



2Tti< 



^1 + iqc"^ + ]_gi^) 



(17) 

See Appendix for a brief derivation. 

The times scales for the orbital decay of plan- 
etesimals or protoplanets under aerodynamical 
and tidal drag are shown in Fig.l (See page 9). 
During the disk evolution and depletion (on a 
time scale ^ 10^^^ yr), aerodynamical drag is 
more effective for small planetesimals with mass 
< 10^'^g, while tidal drag is more important for 
embryos with mass > 10^'^g. The surface den- 
sity of the gas is globally depleted on a time scale 
Tdep ~ 1 — 3Myr, so the magnitude of all time 
scales of the gas drag increases. 

2.2. Protoplanet model 

In this paper, we study the growth of an iso- 
lated protoplanet as the progenitor of Jupiter. 
We analyze the dynamical evolution of its nearby 
residual planetesimals and embryos subject to the 
perturbation of the protoplanet and gas drag. We 
start at the beginning of the quasi-hydrostatic sed- 
imentation stage (S2) of the protoplanet forma- 
tion. The following simplifications are adopted in 
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our simulations: 

i) After its core has acquired an isolation mass, 
the protoplanet is assumed to accrete gas with 
prescribed rates. The gas accretion model is de- 
scribed in §2.2.1. 

ii) The dynamical feedback of planetesimals to the 
protoplanet is neglected because random orbital 
phases usually counteract each other. The accu- 
mulative feed-back perturbations by its co-existing 
nearby embryos on the protoplanet will be studied 
in the future. 

iii) The protoplanet does not have a significant 
radial excursions during the evolution. 

The last approximation is consistent with ne- 
glecting type I migration. Although close-in plan- 
ets may have attained their present-day orbits 
through extensive type II migration (Lin et al. 
1996), these events occur after the proto-gas-giant 
planets have already acquired most of their masses 
and on the viscous evolution time scale of the pro- 
tostellar disks (Lin & Papalozou 1986). Here we 
focus our investigation primarily on the acquisi- 
tion of planetesimals during the formation stages 
of a protoplanet. 

Under the above simplifications, we place a pro- 
toplanet in an orbit with elements Op = 5 AU, 
Bp = 0.01, and ip — 0.005. According to equa- 
tion (fT7|) . the eccentricity and inclination of such 
an isolated protoplanet embedded within the so- 
lar nebula would be damped by the tidal drag 
soon, on a time scale shorter than its growth time 
scale prior to gap formation but longer than the 
synodic period of nearby embryos and planetesi- 
mals (see below). In dynamical equilibrium, the 
recoil motion of the protoplanet is balanced by 
the tidal damping to yield the assumed eccen- 
tricity and inclination, especially during the onset 
of the quasi-hydrostatic sedimentation stage (S2) 
when the mass of the protoplanet is modest. Dur- 
ing the runaway gas-accretion stage (S3), a gas 
gap forms near the protoplanet and Ttidai becomes 
longer than its growth time through gas accretion. 

2.2.1. Inner and outer feeding zones 

Planetesimals grow into protoplanetary em- 
bryos (which significantly perturb the motion of 
their neighboring planetesimals) and cores (which 
accrete gas) through cohesive collisions (Safronov 
1969). The region from where a protoplanet has 



a non-zero probability to accrete planetesimals 
during a single azimuthal passage is referred to 
as its feeding zone. Neglecting the host stars' 
tidal perturbation on their equation of motion 
(but including the Keplerian shear), all planetes- 
imals with semi-major axis a and eccentricity 



e > 5a = |a/ap 



1| are contained in the feed- 



ing zone of the protoplanet centered on its semi- 
major axis Op. However, many synodic periods 
("^syno 2Pk/(35Q), where Pk is the period of 
Kepler motion) may be needed before close en- 
counters can occur. The accretion rate onto the 
protoplanet with a mass Afp and radius Rc is 
(Safronov 1969) 



n2 ( ^GMp 



(18) 



where a = -|- 9Sl/AQ]^t 



Sd/2-ffp, and Sd are the velocity dispersion, scale 
height, spatial and surface density of the planetes- 
imals, respectively. The eccentricities of the plan- 
etesimals are excited by the secular perturbation 
of the protoplanet, and the average excursions of 
the eccentricities and inclinations of the planetes- 
imals per synodic encounter can be expressed as 
(Hasegawa & Nakazawa 1990), 



< Ae >f 



1.9h 



da, < Ai >? 



i.a/i 



(19) 

where h = (Mp/SM^)'^^^ is the scaled HiU's radius 
of the protoplanet. Note that secular perturbation 
does not significantly change the semi-major axes 
of the planetesimals. Over Tgyno the eccentricities 
of planetesimals with mass larger than 10^* g are 
not significantly damped by either aerodynamical 
or tidal gas drag. In principle, nonlinear diffusion 
can lead to further eccentricity growth in a gas- free 
environment. But it proceeds on time scales much 
longer than Tgyno and is effectively suppressed by 
the gas damping effect (Zhou et al. 2007). 

Equation implies e ^ Sa at Sa ^ 1.5ft.. We 
define this location to be the boundary between 
the inner and outer feeding zone. In the inner 
feeding zone, the planetesimals undergo radial ex- 
cursions which cross the protoplanet orbit at each 
azimuthal conjunction. The two-body formula is 
a reasonable approximation to their encounters. 
In the outer feeding zone where Sa > 1.5ft., two- 
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body effects alone cannot lead to close encoun- 
ters. However, when the tidal perturbation of the 
host star is also included in the equation of mo- 
tion, the feeding zone expands to include more 
distant planetesimals. A simple approach to ap- 
proximate the orbits of the planetesimals is to use 
a restricted three-body approximation. Although 
the perturbation of the protoplanet induces peri- 
odic (or chaotic if the planetesimal has relatively 
large energy) variations on the planetesimals' star- 
centric orbital elements a, e, i, the motions of the 
planetesimals are constrained by the Jacobi inte- 
gral, which can be expressed as 

Ej^\{e' + t')-l6l + lh' + 0(h'). (20) 

Planetesimals with positive Jacobi energy re- 
side in the feeding zone (Hayashi et al. 1977), 
i.e. they have finite close encounter probabil- 
ity per Tgyno- In the (a, e) and (a,i) planes, the 
boundary of the feeding zone Ej — is on a 
branch of hyperbolic curves. The half width of 
the feeding zone is Sa ~ 2^/3h for planetesimals 
with e < Sa- In the outer feeding zone where 
5a/h ~ 1.5 — 3.5, e < Sa, the planetesimals can 
engage in close encounters and possible physical 
collision with the protoplanet only within certain 
ranges of the longitude of periapse. Through non- 
linear diffusion, the orbit-orientation angle under- 
goes random walk (Zhou et al. 2007) and over 
many Tsy„o's planetesimals may occasional ven- 
ture into the collision "key holes" of the proto- 
planet. The range of the phase for direct colli- 
sion decreases with Sa and vanishes in the limit 
Sa > 2\^h. Thus, the accretion rate of planetesi- 
mals onto the protoplanet from the outer feeding 
zone with (5a ^ 1.5 — 3.5ft. is significantly reduced 
from the value of Mp in equation (fT5)) . 

However, in a gaseous environment, the accre- 
tion of planetesimals in the entire feeding zone can 
be enhanced by the gas drag effect, especially in 
the outer feeding zone. Although the eccentricity 
excitation of the planetesimals in this location is 
limited by equation (|19p , tidal damping of their ec- 
centricities leads to orbital decay. In i iS.ll we show 
that planetesimals drift in from the outer edge of 
the feeding zone to the proximity of the proto- 
planet where the gravitational perturbation is in- 
tense and close encounters occur more frequently. 
It is such induced orbital decay rate rather than 



the protoplanet accretion rate (Mp) that deter- 
mines the time scale for the clearing of the feeding 
zone. Along the way, these migrating planetesi- 
mals pass through the mean motion resonances of 
the protoplanet and their eccentricities are excited 
to large amplitudes. The subsequent gas drag 
leads to orbital decay. With relative small Taoro 
and Ttidai, small planetesimals and large embryos 
pass through the mean motion resonances. But 
some intermediate-mass planetesimals have rela- 
tively long Tt^a and Ta,a and may be trapped by 
the mean motion resonances (see H3A^ . 

2.2.2. Isolation core mass and gas accretion 
model 

A planetary core temporarily halts its growth 
by accreting planetesimals when it attains an iso- 
lation mass Miso (Lissauer 1987). If all the plan- 
etesimals in the feeding zone can be accreted by 
the protoplanet, the magnitude of Mjso would be 
determined by Sd and the width of the core feed- 
ing zone (Aa): 

Miso = 27ri;dapAa. (21) 

It is useful to scale the magnitude of Sd with 
that of the fiducial minimum solar nebula (Hayashi 
1981) outside the snow line, 

Sd,min - 30(a/lAU)-3/2 g cm-2, (22) 

by a multiplicative factor /d. Numerical simula- 
tions (Kokubo & Ida 2002) indicate that the feed- 
ing zone of a core has a width of = lOop/i 
which is slightly larger than twice 2y/3aph for 
low-eccentricity planetesimals. This minor expan- 
sion is due to the eccentricity excitation associated 
with nonlinear diffusion (rather than linear secu- 
lar perturbation) over the time scale for reaching 
the feeding zone, ~ 3Miso/Mp(Miso). From these 
values, we obtain 

M.o = 2Me/,^/^(^)'^'. (23) 

We place our protoplanet at the present-day loca- 
tion of Jupiter, i.e. at 5AU which is slightly out- 
side the snow line. In a minimum mass nebula, 
this radius is a preferred location for the onset of 
gas giant formation because 1) the isolation mass 
of the embryos is a few and 2) the time scale 
for the buildup of the embryo with isolation mass 
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is comparable to or shorter than the disk depletion 
time scale (Ida & Lin 2004). We take /d = 2 in 
this study, so the protoplanet core in 5 AU has a 
mass of Mpo = 5.67 before the onset of gas 
accretion. 

Cores with isolation masses attract nearby gas 
because their surface escape speed is much larger 
than the sound speed of the disk gas. Neverthe- 
less, heat is released from the contraction of gas 
onto the cores. When the protoplanet mass is low 
(a few Mg), inefficient heat transfer in its enve- 
lope leads to the buildup of a high pressure gra- 
dient to balance its gravity and slow down the 
accretion rate in the quasi-hydrostatic sedimen- 
tation stage (S2). Early numerical models of pro- 
toplanetary structure indicate that the main heat- 
diffusion bottleneck is in the outer radiative region 
(Pollack et al. 1996). The gas accretion rate would 
be greatly enhanced if the grain-dominated opac- 
ity is suppressed (Ikoma et al. 2000, Hubickyj et 
al. 2005). A possible mechanism for major opac- 
ity reduction is dust coagulation and sedimenta- 
tion through the outer radiative region. However, 
the grains may also be replenished by much larger 
colliding planetesimals which disintegrate during 
their passages through the envelope. The clear- 
ing of planetesimals from the feeding zones of the 
accreting cores would greatly reduce the resupply 
rate. 

As a first approximation, we assume the heat 
transfer barrier can be bypassed by the depletion 
of the feeding zone, and the process of gas accre- 
tion is not impeded by the radiative feedback. At 
5AU, when the protoplanet mass is well below Mj, 
its Bondi radius Ri, = GMp/c^ is smaller than its 
Hill's radius Rh — hap and the disk scale height 
H/ap — Cs/V^k = 0.07(ap/5 AU)^/^ in a minimum 
mass nebula. For such a protoplanet, the disk gas 
in the background is homogeneous and the tidal ef- 
fect of its host star can be neglected. For most of 
our calculations, we adopt the conventional Bondi 
formula (Frank et al. 2002) for spherical accretion, 
in which the accretion rate is given as 

Mbd = ^^aM^ = ^aM^ (24) 

where a is a constant of order unity determined by 
the state equation of the gas. For a protoplanet 
at 5 AU, we find ^ ~ 1-2 x lO-^s/gg-iyr-i = 
24/gMQ^yr^^. By integrating equation we 



obtain M = Mo[l - (1 - Mo/Mf){t - to)/rg,.ow]"^ 
for t — to < Tgrow, where Iq is the epoch when 
the accretion begins, Tgrow is the time scale of the 
protoplanet growth to a mass of Aff, 

1/1 1 \ 1 

\Mo Mi J 24/gQ!/zo 

where /^o — Mq/Mq, Mq is the mass of proto- 
planet at the onset of gas accretion, and we assume 
that Mf > > Mq in the approximation of equation 
(I25p . In a minimum mass nebula, a 5.67 pro- 
toplanet at 5 AU has Tgrow ~ 10^ yr. This time 
scale only refers to the growth time scale associ- 
ated with the gas accretion. The emergence of 
such a massive core and the transition from stages 
of embryo-growth (SI) to quasi-hydrostatic sedi- 
mentation (S2) may take longer time. 

The validity of the assumed homogeneous, 
unimpeded, spherical accretion flow onto the 
protoplanet is questionable when its mass be- 
comes comparable to that of Jupiter because 
Rb ~ Rfi ~ H. The effect of differential rota- 
tion in the disk and the tidal torque by the host 
star channel the accretion flow through a proto- 
planetary disk. We are primarily interested in the 
dynamical evolution of residual planetesimals and 
embryos near the outer regions of the feeding zone 
(da > h) which is not strongly perturbed by the 
distribution of gas within R^. As the protoplanet 
attains its asymptotic mass, its tidal interaction 
with the gas disk leads to the formation of a gap 
(see i j2.2.3p . The open of gap reduce the accre- 
tion rate from the unimpeded M ~ 10"'^ Mj yr"-'^ 
in equation (|24|) by several orders of magnitude 
(Dobbs-Dixon et al. 2007). However, this tran- 
sition occurs rapidly and we can approximate it 
with an abrupt termination of its growth. In or- 
der to take into account the uncertainties in the 
boundary conditions, we consider a range of val- 
ues for a and /g so that Tgrow = 10^ — 10^ yr in 
the following calculations. 

In the Bondi model, the growth time scale de- 
creases with the mass of the protoplanet. As we 
show below, this growth pattern can lead to the 
initial quenching of planetesimal bombardment 
and the late-stage capture of residual embryos and 
planetesimals. In order to highlight this behavior, 
we consider a second series of models with an ad 
hoc prescription in which a constant gas accretion 
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rate, 



Mln = 



-(Mf-Mo), 



(26) 



onto the protoplanet is assumed. Figure 2 shows 
the evolution of protoplanetary mass by accreting 
gas on a time scale of Tgrow = 10^ yr according to 
these two models. The initial mass of the isolated 
core is set as 5.67 Mq. Gas accretion is most ef- 
fective around Tgiow = 10^ yr in the Bondi model. 

2.2.3. Gap opening in gaseous disk 

When a protoplanet grows to a sufficiently large 
mass, a gap in the gas disk forms around its orbit 
(Lin & Papaloizou 1979). The formation of the 
gap has two important effects. As we have already 
indicated above, gap formation greatly reduces the 
accretion rate onto the protoplanet and effectively 
terminates its growth. The clearing of the gas also 
reduces the magnitude of the drag on the residual 
planetesimals and embryos (see i 32.ip . Here we 
briefly describe our prescription for the emergence 
of tidally induced gaps in the disk. 

The critical mass of a planet (Mc) over which it 
can open a gap is determined (Lin & Papaloizou 
1993) by a viscous and a thermal condition, 




10" 10" 10" 10'° 10'' 10'' 10" 10'* 10'= 10" 10" 10'" 10" 
mass [g] 

Fig. 1. — Damping time scales in the semi-major axis 
o and eccentricity e due to the aerodynamical (Eq. 
[13]) or tidal drag (Eq.[17]) for different mass of a plan- 
etesimal (or a protoplanet) located at 5 AU. For aero- 
dynamical drag, a density of p = 2 g cm~"^ is assumed. 



Mr 



Mr, 



Mr 



(27) 



(28) 



Mq ' \apj 

where = ac^ / ilk is the kinematic viscosity, and 
a is a parameter (Shakura & Sunyaev 1973). At 
5AU in a minimum mass nebula, a protoplanet 
continues to grow until its mass reaches Mc^t ~ 
300 so that both the thermal and viscous con- 
ditions are satisfied . 

The minimum width of the gap is the Hill's ra- 
dius Rfi of the protoplanet. Numerical simulation 
indicates that the gap extends to 2 — 3Rh. As the 
mass of the protoplanet mass approaches that of 
Jupiter, its unimpeded growth time scale reduces 
to less than 10'^ yr, which is shorter than the vis- 
cous diffusion time scale across all regions of the 
disk wider than Rh- Consequently, gas within a 
gap-feeding zone is rapidly depleted. In the low- 
viscosity limit (where a << 1), we estimate the 
half-width of the gap A on the assumption that all 
the gas in the gap is accreted on to the protoplanet 




2 3 4 5 

log t [yr] 

Fig. 2. — Increase of protoplanet mass in two gas- 
accretion models: Bondi accretion (Eq. [23]) and lin- 
ear accretion (Eq. [21]) on a time scale Tgrow = 10^ 
yr. The initial mass of the planet is 5.67 and the 
final mass is 1 Jupiter mass. 



9 



and contributes to its asymptotic mass. Let us de- 
note the mass increment as 5M . In the absence of 
viscous mass diffusion, 5M = Sn J^^ p^aHda, 
and the half-width of the gap is given by 




From the above minimum-mass nebula parame- 
ters and the present-day Cp — bAU , SM + MpQ = 
1O~^M0, we find A « 0.8 AU, which is approxi- 
mately 2.3Rh and comparable to the value found 
in numerical simulations (cf Bryden et al. 2000). 

The gap-opening in the gas disk near a massive 
protoplanet may change the kinematics of nearby 
planetesimals. At the inner edge of the gas gap, 77 
is enlarged due to the sharp pressure contrast, and 
the orbital decay rate of planetesimals is enhanced. 
However, at the outer edge of the gas gap, the di- 
rection of gas pressure is inwards and the pressure 
gradient is positive, resulting in > accord- 
ing to equation ([J), and rj is negative by equation 
(HI) (Lin & Papaloizou 1993, Bryden et al. 2000). 
Thus the tail wind of the gas induces the planetes- 
imals to migrate outward. The severe depletion of 
gas in the gap region also reduces the drag at that 
location. However, these feedback effects are neg- 
ligible until the protoplanet has opened a gap and 
attained most of its asymptotic mass. After the 
gap formation, the sign of a becomes positive (see 
Eq. [13| ) for low-eccentricity planetesimals (with 

< |ry|) and they move beyond the gap, whereas 
the high-eccentricity planetesimals continue to mi- 
grate inward. 

The modification of Vg near the outer edge of 
the gas gap also locally changes the sign of the 
tidal torque on the protoplanetary embryos (see 
Eq. [15]). In principle, this process induces the 
embryos to undergo outward migration, analogous 
to the type I inward migration. But similar to 
the aerodynamical drag, the damping of the em- 
bryo eccentricities dominates the evolution of their 
semi- major axes. 

As the depletion of gas in the gap reduces the 
local density of the gas disk, we take f = in the 
corresponding aerodynamical and tidal drag for- 
mulas (Eqs. [1] and [15]) inside the gap in our 
simulations. We assume that gas in the gap is de- 
pleted on a time scale r = Ar'^/iy ~ H^/i^, which 
gives r « 2a~^ yr. Thus for a = 10~^, we have 
r « 2000 yr. However, to simplify the problem, 



we neglect the migration effect due to gap forma- 
tion, which is very interesting and will be studied 
elsewhere. 

2.3. Planetesimal model 

In this paper, we study the evolution of plan- 
etesimals with initial semi-major axes a G [3.2, 8.5] 
AU, i.e., a G [cip — 5ap/if,ap -I- lOop/if] with Op = 
5AU, h{ = (Atf/3)i/3 and /Xf = lO'^. We choose 
a larger region for the outer part, because due to 
gas drag and eccentricity damping, planetesimals 
on both side of the protoplanet suffer orbital de- 
cay. Exterior to the protoplanet, planetesimals 
migrate toward its feeding zone. 

Numerical simulations indicate that during the 
oligarchic growth, the feeding zone of the most 
massive embryos may indeed contain less massive 
planetesimals and embryos (Kokubo & Ida 2002). 
In order to examine the mass-dependent collision 
probability of planetesimals with a growing proto- 
planet, we carry out several simulations, each with 
a population of uniform mass planetesimals. The 
individual planetesimal mass under investigation 
ranges from 10^^ — lO^^g. The planetesimals are 
treated as test particles in the simulation (i.e. they 
do not impose any gravitational perturbation on 
each other or on the protoplanet) except when we 
evaluate the magnitude of the gas drag. The low- 
mass range corresponds to km-size objects, which 
can withstand the aerodynamical drag in a mini- 
mum mass nebula with Ta^a > Tdep and be retained 
by the disk. We did not simulate the interaction 
between a protoplanet and embryos more massive 
than 10^''g because it would not be adequately ap- 
proximated by the restricted three-body approach. 
A more comprehensive treatment of the interac- 
tion between a population of comparable mass em- 
bryos will be presented elsewhere (see Zhou et al. 
2007). 

We choose three sets of models for detail anal- 
ysis. In Models 1, 2, and 3, we set the mass 
of each planetesimal to be m = 10^^ g (Model 
1), 10^4 g (Model 2), and 10^^ g(Model 3). The 
planetesimals in Model 2 correspond to the tran- 
sitional objects as planetesimals evolve into oli- 
garchies, which perturb the velocity dispersion of 
their neighbors. The embryos in Model 3 corre- 
spond to the isolation mass in a minimum mass 
nebula interior to the snow line. These three rep- 
resentative models are also chosen to illustrate the 
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relative importance of the gas drag effects. The 
dominant physical process for eccentricity damp- 
ing is aerodynamical drag for the low-mass plan- 
etesimals in Model 1 and tidal drag for the em- 
bryos in Model 3. The eccentricity damping time 
scale is the longest for the intermediate-mass oli- 
garchies in Model 2. 

In order to build up sufficient samples for a 
statistical analysis, we use 1000 planetesimals for 
each simulation. We normalize these models with 
= /d5]d,min, SO the total mass of planetesimals 
in the region [ain, flout] is: 



Mtot = XTj""' /d Sd,min27rada 



>l/2 



.1/2 



lAU^ 



.1 AU/ 



(30) 

Thus the total mass of the planetesimals in 
[3.2, 8.5] AU with /d = 2 is 31.5 M®. 

Due to their interaction with each other and 
the turbulent gas, the velocity dispersion of plan- 
etesimals, = V — Vk, is expected to have 
a Gaussian distribution in Cartesian coordinates. 
The corresponding eccentricities and inclinations 
of the test particles follow a R ayleigh distribution 
( Greenzweig fc Lissauer 1992[ l. The initial aver- 
age eccentricities and inclinations of the planetesi- 
mals are taken as 0.0007 and 0.00035, respectively. 

3. Formation of planetesimal gap during 
gas accretion 

For our numerical simulations, we use the Reg- 
ularized Mixed Variables Symplectic (RMVS3) 
method in the SWIFT packet (Levison & Duncan 
1994). This algorithm is well suited for the com- 
putation of close encounters of planetesimals with 
protoplanets (but not the inter-particle encoun- 
ters). The time step of the integration is set as 0.05 
yr, i.e., ~ 1/200 of the period of the protoplanet 
orbit. All planetesimals which venture within the 
physical radius of the protoplanet are assumed be 
accreted by it. The density of planetesimals is 
taken as ptp = 2 g cm^'^ in equation (fT2|) . The ra- 
dius of the protoplanet is given by equation pT] ) 
with two limiting values of pp = 1 g cm~^ and 
0.125 g cm^''. The quantity pp — 1 g cm~^ cor- 
responds to that of the core with radius Re- In 
reality, the effective capture cross section of a pro- 
toplanet is determined by the planetesimal masses, 
the relative speed as well as the protoplanet mass 



(which determines the density distribution in its 
gaseous envelope). During the initial gas accretion 
stage of a protoplanet, its envelope has a relatively 
small mass and does not significantly modify the 
capture cross section. After the onset of efficient 
gas accretion, the density in the outer regions of 
the envelope is also low as it undergoes dynami- 
cal collapse. Gas accumulates near the core as its 
inward motion is halted. The radius of this loca- 
tion is determined by the efficiency of radiation 
transfer but is a few times larger than Rc. After 
the protoplanet has attained its asymptotic mass 
and gas accretion onto it is quenched, the radius 
of a typical gas giant quickly reduces to twice that 
of Jupiter. In order to take these possibilities into 
account and still keep our investigation within rel- 
atively general and simple bounds, we consider a 
second set of capture criterion for intruding plan- 
etesimals by doubling the protoplanet 's effective 
radius (which corresponds to assuming a homoge- 
neous density pp — 0.125 g cm~'^ for the proto- 
planet). 

There are also planetesimals scattered to solar 
distances larger than 100 AU or smaller than 1 
AU. These particles are respectively classified as 
ejectors to the outer solar system or Sun crashers. 
The planet mergers, ejectors, and Sun crashers are 
registered and removed from the subsequent evo- 
lution. 

The eccentricity of the planetesimal orbits are 
excited by the protoplanet and damped by gas 
drag. Eccentricity damping also lead to a decline 
in the semi- major axes (see §2.ip . The energy of 
the planetesimal orbits may also be modified by 
resonant interaction and close encounters with the 
protoplanet. Planetesimals with a < 1 AU are 
not able to collide with the protoplanet in their 
subsequent evolutions, therefore an inner bound- 
ary of 1 AU is adequate for the objectives of the 
present investigation. In all the simulations pre- 
sented here, we assume the gas disk has an initial 
surface density of the minimum solar nebula (see 
Eq. [9]). After the formation of the protoplanet 
core, a uniform exponential depletion of gas on a 
time scale of 10^ yr is also assumed. In most but 
not all models, this time scale is larger than the 
magnitude of Tgrow The solid disk has an initial 
surface density of twice the minimum solar nebula 
(see Eq. [22]). 
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3.1. Opening of planetesimal gap 

An important phenomenon in the evolution 
of planetesimals in a gaseous environment is the 
opening of a planetesimal gap around the pro- 
toplanet. The opening of a planetesimal gap in 
a gas-free environment was discussed by Rafikov 
(2001). This situation is analogous to planetary 
rings being shepherded by a satellite. Due to the 
differential Keplerian motion, inter-particle en- 
counters lead to angular momentum transfer and 
the diffusion of the ring. But in the proximity 
of the protoplanet, tidal perturbation from the 
protoplanct tends to expel the planetesimals away 
from it (Goldreich & Tremaine 1978, 1980, Lin & 
Papaloizou 1979). When the tidal torque exceeds 
the rate of inter-particle angular momentum ex- 
change, a gap centered on the protoplanet forms. 

During the growth of the embryos and the for- 
mation of the core, the residual planetesimals at- 
tain a MRN size distribution, i.e. most of the plan- 
etesimals' surface area and mass are contributed 
by the small and large planetesimals respectively. 
During the stages of quasi-hydrostatic sedimen- 
tation (S2) and runaway gas-accretion (S3), the 
maximum embryo size has already increased, so 
that the collision frequency of the smaller resid- 
ual planetesimals becomes small compared with 
that of their synodic encounters with the proto- 
planet. In this low-collision- frequency limit, the 
secular perturbation of protoplanet on the resid- 
ual planetesimals needs to be taken into account 
(Franklin et al. 1980, see pXT|l . 

Although the energy dissipation rate associated 
with planetesimal collisions is reduced with the 
collision frequency, the excited planetesimals also 
experience eccentricity damping from the disk gas. 
Since Ta,e and tj e in equations (fT3|) and (fT7|l are 
much larger than Tgyno (over which span the eccen- 
tricities of planetesimals are excited to < Ae >), 
the gas drag does not directly reduce the plan- 
etesimal eccentricities from that in equation (jl9p . 
However, this process is accompanied by a slow 
radially inward drift on time scales Ta,a ^ e.~'^Ta,e 
and Tt^a ^ e~^T(_e (see Eqs. [13] and [17]) . This 
drift speed is faster for planetesimals inside than 
those outside the feeding zone of the protoplanet 
because they are excited to larger < Ae >. 

The inward drift causes planetesimals interior 
to the protoplanet orbit to drift away from the 



feeding zone of the protoplanet. But it also brings 
the planetesimals exterior to the protoplanet orbit 
into its feeding zone, where their orbital responses 
become chaotic, resulting in capture or close en- 
counters. Planetesimals on both sides of the pro- 
toplanet orbit are evacuated and a gap eventually 
appears. Beyond the feeding zone, the planetesi- 
mal orbits evolve slowly because the protoplanet 
perturbation is much weaker so that the magni- 
tude of < Ae > is much smaller. Some external 
planetesimals are trapped onto mean motion res- 
onances of the protoplanet, where their eccentric- 
ities are increased. Since the gas gap is confined 
to within the feeding zone, the residual disk gas 
damps the planetesimal eccentricities to modest 
equilibrium values. 

Figure 3 illustrates the distribution of low-mass 
(m = 10^^ g) planetesimals (Model 1) in the (a, e) 
plane at different epochs. In this model, we adopt 
the Bondi gas-accretion prescription and set the 
time scale of Tgrow = 10^ yr. Kt t = 10^ yr, 
Afp ^ 6M0, Tgyno ^ 10"^ yr, and the eccentrici- 
ties of the planetesimals within the current feed- 
ing zone are excited to values < Ae >~ 0.01 — 0.1 
(Fig. 3a). After t = 10^ yr, the protoplanet ac- 
quires its full asymptotic mass. There are many 
planetesimals inside the final feeding zone of the 
protoplanet. A planetesimal gap begins to form, 
albeit at a slower rate than the expansion of the 
hypothetical feeding zone (Fig. 3b). At the outer 
regions of the feeding zone, the eccentricities of 
many planetesimals are greatly excited. The V 
shape of the planetesimals' (a, e) distribution in- 
dicates that many are scattered analogous to the 
KBO's. After t > 10^ yr, planetesimals in the 
feeding zone are completely cleared, as those in- 
terior to the protoplanet orbit drift inward while 
those exterior to it are scattered outward (Fig. 3c, 
3d) . The analysis for the clearing time scale of the 
planetesimal gap is presented in §3.21 

Since the opening of the planetesimal gap is the 
result of gas drag and the protoplanet perturba- 
tion, the time scale for their eccentricity damping 
and therefore gap-clearing depends on their mass 
(Ida & Lin 1996). Figure 4 displays the configura- 
tions of survived planetesimals in the (a, e) plane 
for planetesimals with mass m — 10^^ g (Model 
2) and protoplanetary embryos with m = lO^'' g 
(Model 3), respectively. 

During the buildup of the protoplanet asymp- 
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Fig. 3. — Configurations of survived planetesimals 
witli mass m = 10^^ g in tfie {a, e) plane during evo- 
lution. The protoplanet accretes gas according to the 
Bondi model on a time scale Tgrow = 10^ yr. The dot- 
ted and solid lines are the boundaries of the feeding 
zone {E = 0) at time t = 10* yr and t = 10^ yr, respec- 
tively. The gray dot at a ~ 5 AU shows the position of 
the protoplanet. Some locations of mean motion reso- 
nances between planetesimals and the protoplanet are 
shown in plot (d). 
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Fig. 4. — Configurations of survived planetesimals 
with mass m = 10^'' g and m = 10^^ g in the (a, e) 
plane during evolution. The protoplanet accretes gas 
by the Bondi model on a time scale Tgrow = 10^ yr. 
The dotted lines and solid lines are the boundaries 
of the feeding zone {E — 0) at time t — 10* yr and 
t = 10^ yr, respectively. The cases of t = 10* yr are 
similar to Fig. 3. The gray dot at a = 5 AU shows the 
position of the protoplanet. 



totic mass (at t = 10^ yr in Figs. 4a and 4c), resid- 
ual planetesimals are found inside the feeding zone 
in all models. A closer inspection indicates that 
the opening of the planetesimal gap is less efficient 
in Model 2 than in Models 1 and 3. This minor dif- 
ference supports the conjecture that the clearing of 
the gap is due to the combined action of the pro- 
toplanet perturbation and gas drag, because the 
e-damping and a-decay rates of the intermediate- 
mass planetesimals (represented by Model 2) are 
smaller than those of the low-mass planetesimals 
and high-mass embryos (see Fig.l for the drift time 
scale). On time scales much longer than Tgiow, the 
intermediate-mass continues to occupy the edge 
regions of the feeding zone. Many planetesimals 
are also trapped in the outer mean motion res- 
onances with increased eccentricity (Fig. 4b). In 
contrast, the drift speed of the embryos (m = 10^^ 
g) is faster than those in Models 1 and 2, and all 
the embryos in the feeding zone are cleared out 
after t = 10^ yr (Fig.4d). Nevertheless, the width 
of the planetesimal gap is limited to ~ 2\/3hap. 

The opening of a planetesimal gap near the 
protoplanet is also shown in Fig. 5 with the 
intermediate-mass (10^'^ g) planetesimals. At 
t = Tgrow = 10^ yr, the radial distribution of 
the planetesimals show a diffusion profile around 
the corotation radius of the protoplanet. Some 
survival planetesimals near the protoplanet are 
caught onto horseshoe or tadpole orbits. Figure 
6 plots two examples in the tadpole orbits librat- 
ing around the L4 and L5 points, respectively. In 
Fig. 5, there is a slight enhancement of surface 
density beyond the edge of the feeding zone due 
to the resonant trapping of inward-drifting plan- 
etesimals, which will be discussed in §4.3. This 
accumulation of planetesimals increases the local 
surface density and Mjso, promotes the growth 
rate of protoplanetary cores and the emergence of 
secondary proto-gas-giant planets. 

3.2. Planetesimal gap clearing time scale 

There are several relevant time scales to be con- 
sidered. In t j2.2.11 we showed that, under the pro- 
toplanet perturbation, planetesimals with overlap- 
ping 2-body (neglecting the stellar tide) orbits ex- 
tend to the boundary of the inner feeding zone 
where 6a < 1.5/i. The time scale of depletion of 
the total planetesimal population in this region, 
Mp.f — AirJ^dSaoi, can be derived in view of equa- 
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Fig. 5. — Semi-major axes distribution of survived 
planetesimals at four epochs during the mass growth 
of the protoplanet. Each bar corresponds to 0.1 AU. 
The mass of individual planetesimal is 10^^ g. The 
protoplanet accretes gas by the Bondi model on a time 
scale Tgrow = 10^ yr. A^o is the total number of plan- 
etesimals survived at that time. The solid curve cor- 
responds to the initial density profile. The opening of 
a planetesimal gap leads to a slight enhancement of 
surface density near the boundary of the feeding zone. 
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Fig. 6. — Tadpole orbits of two surviving planetesi- 
mals in model Fig.Sd. The coordinate is set as origin 
at the mass center of the sun and protoplanet (with 
mass growing), and corotating with the protoplanet at 
~ 5AU. The planetesimals' asymptotic eccentricity is 
e ~ 0.001. The dot around (5, 0) marks the orbit of the 
protoplanet. The planetesimals have semi-major axis 
a = 4.875 and a — 5.038 at t — 10^ yr, respectively. 



tion (fT8|) with e ^< Se > from equation p9|) . 




where Pk is the Keplerian period. So Tp^f < 10^ 
yr. 

This rapid depletion time scale is only applica- 
ble to the planetesimals in the inner feeding zone 
where Sa < l.bh. Planetesimals in the outer feed- 
ing zone with Sa/h ^ 1.5 — 3.5 only occasionally 
venture into the Roche lobe of the protoplanet. 
Based on the discussions in the previous section, 
we now derive the time scale for planetesimals to 
migrate across the boundary between the inner 
and outer feeding zone. Since Tp^f is relatively 
small, the duration of the migration across the 
outer feeding zone essentially corresponds to the 
clearing time scale of the planetesimal gap. 

Suppose the protoplanet has a mass /z = 
Afp/M* at time t, so its instantaneous normal- 
ized Hill radius is h = In terms of 
the protoplanet 's asymptotic (at t > Tgrow) nor- 
malized Hill's radius /if, the scaled distance of a 
planetesimal is defined as 

where 5a = \a/ap— 1|. According to equation (|13|). 
the speed of the inward drift for a planetesimal 
with mass m < 10^'^ g is given as 

aacro ~ -1.7^6^, (33) 
Tacro 

where we suppose » 77, because the eccentric- 
ity of planetesimals inside the feeding zone could 
be excited up to ~ 0.1 (e.g., Figs.3-4). In ^Z2l[ 
we evaluate the average excursions of < Ae > 
per encounter. (As small initial inclination is ex- 
pected, we neglect the modulation in Ai.) Since 
Ta,e >> Tsyno, wc Can replace e in equation ([551) 
by < Ae > in equation to obtain 

550 1 

(6f)aero=- Tg T7 ' (34) 

'aero f^f 

where the dependence on h{ is introduced for the 
purpose of normalization. With the growth of 
the protoplanet mass, h also increases with time. 
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However, in order to simplify the problem, we first 
assume h as a constant hf (an approximation to 
be justified a posteri) so that 

(35) 

where /3i = 3.2. The limiting value (2^/3) for 
b{ corresponds to the width of the entire feeding 
zone. Planetesimals drift from the protoplanet 
orbit to the inner boundary of the asymptotic 
feeding zone and from the outer boundary of the 
asymptotic feeding zone to the protoplanet orbit 
on a time scale of 

The above expression is obtained by equating the 
left hand size of equation ([55)1 with 2y/3h/hf. 

According to the above equation, gap formation 
for the low-mass planetesimals proceeds on a time 
scale Taero/h'^- Substituting Tacro from equation 
([H)) and assuming a constant protoplanet mass, 
we find Topen ~ 3(M,7/Mp)2/3 Myr for Model 1. 
We carry out two comparisons, Models la and 3a, 
in which the mass of the protoplanet is fixed to 
that of Jupiter so that h = hf. We compare, in 
Figs. 7a and 7c, the results of the numerical sim- 
ulations with that in equation ([35|l . Numerically, 
we determine hf from the distribution plots such 
as Fig. 5 at some typical epoch. The magnitude of 
hf is defined to be the maximum half- width of the 
gap within which only planetesimals on horseshoe 
or tadpole orbits (around 5AU) survived. 

The qualitative agreement between the func- 
tional dependence of hf on h in the numerical re- 
sults and the expression in equation (|35p supports 
the interpretation that the clearing of the gap is 
regulated by the orbital decay of the planetesimals 
in the feeding zone. (This agreement is particu- 
larly good during the expansion of the gap through 
the initial outer feeding zone where the migration 
scenario is most applicable.) However, in compar- 
ison with the expression in equation ([55]) , a factor 
of 2 — 3 for j3i is needed to fit the numerical sim- 
ulations. This difference in the magnitude of j3i 
is caused by an underestimate in the analytical 
expression for < Ae > which did not take into ac- 
count the cumulative eccentricity modulation dur- 
ing each secular cycle. From equation (I36p , the 



expansion of the gap is stalled with an asymptotic 
width ~ 2y/2>h at Tacro in the numerical sim- 
ulations. A similar analysis also applies to the 
massive embryos for which the tidal drag is more 
appropriate. 

We now return to the more realistic models 
in which the planctcsimal gap formation proceeds 
during the growth of the protoplanet. Despite the 
increases of the protoplanet mass, the above ap- 
proximations would essentially be valid if Tgiow > 
Topen- In this limit, planetesimal gaps with the in- 
stantaneous feeding zone width (2^/3/l) fo rm and 
expand with the mass of the protoplanet. Based 
on the assumption that most planetesimals in the 
feeding zone can collide with the protoplanet core 
and are the main contributors to the initial growth 
of the protoplanet. Pollack et al. (1996) derived 
their bombardment rate onto the core from the ex- 
pansion rate of its feeding zone. In that model, the 
suppression of the gas accretion rate due to the en- 
ergy dissipation of planetesimal accretion has been 
taken into account. 

However, in the limit of modest Mp (a few M®), 
both the protoplanet 's h and the planetesimal ec- 
centricities are very small so that Topcn > Tgto-w 
even for a protracted protoplanetary growth. In 
this case, the feeding zone expands faster than 
they can be cleared out (especially in Model 2 
in Fig. 4). Both the protoplanet mass and feed- 
ing zone width attain their asymptotic values on 
a time scale Tgi-ow, after which gap formation and 
clearing of the planetesimals in the feeding zone 
proceed on a time scale Topon- This gradual mass 
ramp up significantly delays the clearing of the 
gap. During the advanced stage, the assumption 
of constant Mp is again satisfied and the values of 
h can be replaced by hf in the above equations. 
In the next section, we show that not all the plan- 
etesimals in the feeding zone collide with the pro- 
toplanet core and their collision rate may be sub- 
stantially lower than that estimated by Pollack et 
al. (1996). This effect can reduce the barrier for 
the gas accretion rate onto the protoplanet enve- 
lope. 

Using the same approach, we deduce the width 
and time scale associated with the gap-opening of 
intermediate and high-mass (> 10^"^ g) embryos. 
In this case, the gravitational tidal drag provides 
the dominant eccentricity damping effect. From 
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equation p7|) . we find 



5 a 



atidal 



and 



Ttidal 



(38) 

where /32 = 2.7. The time scale to open a gap is 
given as, 



(37) 



opon^tidal 



5 tidal 
^ /32 ^ h ' 



(39) 



Figure 7 shows the evolution of the width of the 
planetesimal gap determined from our numerical 
simulations of Models 1 and 3. In Figs. 7c- 7d, the 
protoplanet accretes gas according to the Bondi 
accretion prescription with Tgiow — 10*^ yr, as in 
the cases of Figs. 3-5. The sohd curves are the- 
oretical predictions given by equations ()35p and 
with different coefficients. The functional 
forms again are in general agreement, though there 
is a factor of 3 discrepancy in the coefhcient of (32 ■ 
Outside the boundary of the feeding zone, the 
rate of the eccentricity excitation by the pro- 
toplanet is small. According to equations 
and (jl7p . the migration speed becomes corre- 
spondingly smaller by several orders of magnitude. 
Though the migration speed outside the gap can 
be obtained following similar lines of reasoning, 
such an analysis is not crucial for this study and 
we will not discuss it further. 

3.3. Dependence on gas accretion model 
and time scale 

In the above simulations, the time scale of gas 
accretion onto the protoplanet, Tgrow: is set to be 
10^ yr. In order to determine the dependence 
of the collision efficiency of planetesimals on the 
growing time scale of the protoplanet, we adopt a 
range of Tgiow, from 10'^ to 10^ yr, with both the 
Bondi ((24|) and linear ((26|) prescriptions for gas 
accretion. Figure 8 illustrates the collision (Pcoi) 
and escape (Pose) probabilities of intermediate- 
mass (10^'^ g) planetesimals (in Model 2) during 
the growth of the protoplanet. Since the evalu- 
ation of the protoplanet radius is based on the 
present density of gas giants, Pcoi should be re- 
garded as a lower limit. Although protoplanets 



have extended envelopes, most of their mass re- 
sides in the core and the density in the envelope 
decreases rapidly with radius. Small and modest- 
size planetesimals may be captured by the proto- 
planet when they pass through its envelope. But 
the large embryos can only merge with the proto- 
planet through direct collisions with the core. 

Throughout the evolution, Topon > igrow, 
so that the expansion of the feeding zone out- 
paces the clearing process (Figs. 3-4). Conse- 
quently, major epochs for planetesimal collisions 
with the protoplanet occur around t = Tgrow 
in Bondi gas-accretion model (Fig. 8). The du- 
ration of the epoch of intense bombardment is 



)±0.25 



for Tgrow > 10 yr. In con- 



trast, major collision events would occur much 
earlier if the protoplanet mass increases according 
to the hypothetical linear gas-accretion prescrip- 
tion. Most of these collisions occur &i t < Tgiow 
The dichotomy between these two types of accre- 
tion arises because the Bondi prescription leads 
to a runaway process, in which most of the pro- 
toplanet mass is attained only at the very end 
when t = Tgrow However, the hypothetical linear 
accretion introduces a much earlier ramp up in 
the protoplanets mass. Consequently its physical 
radius, gravitational perturbation, and width of 
feeding zone also grow quickly, inducing an earlier 
phase of intensified collisions. 

Also note that both Pose and Pool are normal- 
ized to the initial planetesimal population in the 
entire computational domain which covers twice 
the width of the asymptotic feeding zone of the 
protoplanet. The total cumulative statistics sug- 
gest that, the fraction of the original planetesimal 
population in the feeding zone which collides with 
the protoplanet is comparable to that scattered 
to the outer disk regions. With both prescrip- 
tions, the ejections of planetesimals occur around 
t > Tgrow For a given density, the radius and 
surface escape speed Vi^sp of the protoplanet are 

1/3 

proportional to Mp . During the initial growth 
stages of the protoplanet, its I4sp is small com- 
pared with the local Keplerian velocity. Scatter- 
ing from grazing encounters excite the planetesi- 
mal eccentricities rather than eject them. At an 
advanced growth stage of the protoplanet (when its 
Mp ~ Mj), however, scattering with impact pa- 
rameter larger than a few planetary radii can lead 
to large-angle deflections and the escape of plan- 
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Fig. 7. — Evolution of tlie widtii of tiie planetesimal 
gap. Tiie protoplanet lias constant mass — 10~^ in 
panels (a) and (b), and grows from fi = 1.7 x 10~^ to 
H = 10"'' in panels (c) and (d) through gas-accretion 
with a Bondi gas-accretion model of Tgrow ~ 10^ . The 
solid curves are theoretical predictions given by equa- 
tions ((35]) and JSH]) with coefficients: (a) /3i = 7.4, (b) 
132 = 8.0, (c) f3i = 9.6, (d) f32 = 8.0. The squares 
and triangles denote the inner and outer boundaries, 
respectively. The dashed curves in panels (c) and 
(d) show the width of the feeding zone defined by 
a — ap\/{aph) at time t. 
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Fig. 8. — Distributions of the planetesimal collision 
probability Pcoi (a-c) and escape probability Pesc (d-f) 
as a function of the evolution time. The time scales for 
gas accretion are Tgrow ~ 10^, 10^, 10® yr, respectively. 
In the Bondi gas-accretion model, the planetesimals 
collide with the protoplanet mainly around the late 
stage of gas accretion, which is quite different from 
that in the linear gas-accretion model. 



etesimals (Lin & Ida 1997). The ratio of Pcsc/fcoi 
would increase with the semi-major axis (op) and 
effective density (pp) of the protoplanet (Ida & Lin 
2004). The latter quantity is likely to increase 
after the gas accretion onto the protoplanet enve- 
lope is quenched by its tidal barrier (Dobbs-Dixon 
et al. 2007). 

We now scale our model with the minimum 
mass nebula model. By adopting the same ra- 
dial dependence but twice the magnitude in sur- 
face density as that in the minimum mass nebula, 
the total solid mass in the region a E [3.2, 8.5] AU 
is 31.5 Af© in our models(§2.3). We deduce the 
rate of accretion by multiplying this total solid 
mass with the planetesimal accretion probability. 
Figure 9a shows the solid mass (Mcoi) that collides 
with the protoplanet as a function of its growth 
time scale (rgj-ow)- With Bondi accretion, the 
magnitude of M^oi attains a maximum value with 
T G [10^10*^] yr (Fig. 9b), which is 6 - 7 M©. 
(Once again, these values are applied to compact 
protoplanets and should be regarded as lower lim- 
its.) With the same Tgrow, the collided solid mass 
in the linear accretion prescription is 1 ~ 2 M© 
less because the more rapid initial growth of the 
protoplanet leads to early excitation and evacua- 
tion of the feeding zone. This rapid initial ramp up 
of the protoplanet mass causes many planetesimals 
to drift inward due to gas drag or be ejected be- 
fore the feeding zone attains its asymptotic width 
(Fig. 9c). In the limit of large Tgrow, a protoplanet 
with Bondi gas-accretion model can accrete plan- 
etesimals more effectively. 

Figure 9d shows the evolution of solid (plan- 
etesimals) and gas (including dust grains) com- 
ponent accreted by the protoplanet. The proto- 
planet has more solid mass (including the initial 
core which is 6M0) in the beginning of Bondi 
accretion due to the inefficient gas accretion. But 
after Mp = Mgod + ^^gas > 7Af©, the gas accretion 
rate far exceeds that of the planetesimal accretion 
rate because Topcn > Tgrow ■ While the accreted gas 
may carry small dust particles with it, the suppres- 
sion of initial planetesimal accretion reduces the 
feedback effect which limits the gas accretion rate. 
This late addition of planetesimals can lead to the 
enrichment of protoplanet envelope to super-solar 
metallicity. 

We now examine the dependence of planetes- 
imal accretion efficiency on their mass m. This 
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Fig. 9. — Dependence of collided planctcsimal mass 
on Tgrow and gas accretion model, (a) Evolution of 
the planetesimal mass colliding with the protoplanet 
during its gas accretion, (b) Variation of collided solid 
mass with the time scale of gas accretion Tgrow in the 
two different models, (c) A comparison of the sur- 
vived planetesimals with mass m = 10^^ g in Bondi 
and linear gas-accretion models during the evolution, 
(d) Evolution of solid and gaseous compositions of the 
protoplanet. In this paper, wc suppose the solid and 
gaseous disks have a surface density two times and one 
time the minimum solar nebula, respectively. 
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Fig. 10. — Dependence of collision rate on planetes- 
imal mass (m). (a) Evolution of the total solid mass 
collided onto the protoplanet with individual planetes- 
imal mass ranging from lO^'^ — lO^'^ g in the Bondi gas- 
accretion model of Tgrow = 10^ yr. (b) Dependence 
of the collided planetesimals' mass on their individual 
masses. 



dependence arises because the eccentricity damp- 
ing time scale is a function of m (see Fig.l). Fig- 
ure 10 shows the evolution of collided solid mass 
with individual planetesimal mass ranging from 
lO^'^ ~ lO^'^ g. The maximum collision rate occurs 
with individual planetesimal mass m G io20~26 
g, which corresponds to planetesimals with ra- 
dius of 50 ~ 5000 km. These planetesimals also 
have the largest Taero and rtidai- The planet with 
smaller density (p = 0.125 g cm~'^) has a bigger 
collision rate. In Fig. 11, we plot the probabili- 
ties of the planetesimals survived, collided, ejected 
or crashed after 2 x 10^ yr in a Bondi accretion 
model with tgrow = 10^ yr. We find that at least 
2/3 of planetesimals have survived at an epoch 
^ lOtgrow In this model, the initial width of the 
planetesimal disk is Ibuphf out of which a plan- 
etesimal gap with a width ~ Ta^hf is evacuated. 
Although the final distribution of the planetesimal 
disk is more extended (see Fig. 5), most of the sur- 
viving planetesimals (> 2/3) are located within a 
range of nearly Saphf from the protoplanet. Thus, 
the density outside the feeding zone is slightly en- 
hanced on average after the protoplanet obtains 
its asymptotic mass. 

3.4. Resonant trapping 

In both Figures 3 and 4, the accumulation of 

planetesimals near the outer mean motion reso- 
nances of protoplanet is noticeable, particularly 
for the intermediate-mass Model 2. The first-order 
mean motion resonances are located both interior 
and exterior to Op where the period ratio can be 
expressed either as {p+ 1) : p or p : {p+ 1). The 
corresponding ratio of semi-major axes with Op is 
a = (1 + 1/p)^^^'^ or (1 -|- l/p)"^/^ respectively. In 
each case, the protoplanet mass grows from 5.6M0 
to 300 Af0. With its asymptotic mass, the 2:1 
and 3:2 resonances of the protoplanet are outside 
its feeding zone, whereas mean motion resonances 
up to 10:9 are beyond 2^/^h at the onset of the 
simulations. 

Using a simple pendulum model (§8, Murray & 
Dermott 1999), the width and libration time scale 
associated with these mean motion resonances can 
be derived as 

^^^^^^mA^^ap, (40) 

Tres = (3p2eresiir(a)/x)"^^^ flc, (41) 
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where Pk is the Keplerian period, and the magni- 
tude of Rr{a) = a[p+l + ^aD]b^^^ is an increas- 
ing function of p. The minimum value of eccen- 
tricity in the resonance is < Ae > (see Eq. [19]). 
During the resonant passage, an adiabatic invari- 
ant constrains the eccentricity change to be 



Ae 



2 

res 



pap 



(42) 



In the hmit that Aeros due to resonant perturba- 
tion is larger than < Ae > due to secular pertur- 
bation, it can be substituted by ejcs so that 



A, 



3p2 



2/3 



_P_ 

„l/3' 



(43) 



(44) 



y 3 J p^ 

Tres = 12-1/3 ipRr{a)^i)-^/^ Pk. (45) 

Substituting Cres into equations and (fT7)) . 
we find that the damping of the resonantly ex- 
cited eccentricity leads to a characteristic migra- 
tion time across Aios, which is 



^(^rcs) 



/ 16pRr{a)p.'\^ 1/'^ 
V 3~ 



(46) 



for small planetesimals. For the more massive em- 
bryos, 

^rcs 8 

Tx = -. 7 T = -pTtide, (47) 

fltidc(,6rosJ O 

which is independent of the protoplanet mass. 

Planetesimals are trapped in resonances when 
Tx > Trcs- As a planetesimal approaches a proto- 
planet, it encounters resonances with increasing p 
and Rripi), so the magnitude of Tres decrease with 
little change in the magnitude of eres- In principle, 
it should be easier for the strong resonance close to 
the protoplanet to capture the planetesimals be- 
cause Tx becomes larger than Tres for sufhciently 
large p's. However, for relatively large p's, the 
normalized distance separating the p : {p -\-V) and 
(p — 1) : p mean motion resonances, 



2 



(48) 



is a decreasing function of p^ but independent of 
the protoplanet mass. In contrast, equation (|44p 



indicates that the magnitude of Aj-es/op increases 
with /i = Mp/M* and decreases with p^/^. When 
the protoplanet attains 



(49) 



the width of its p : (p -I- 1) resonance (i.e. Aj-es) 
becomes larger than the separation between it and 
the (p— 1) : p resonance (i.e. Ap^p_i). Overlapping 
resonances generally lead to dynamical instabili- 
ties which excite the eccentricities of the trapped 
planetesimals and modify their semi-major axes 
(Murray & Dermott 1999). 

The expression in equation indicates that 
the critical mass for the overlapping resonances 
is a decreasing function of p. During the growth 
of the protoplanet, planetesimals located in the 
initial inner feeding zone become destabilized 
and collide with the protoplanet first. But the 
planetesimals captured onto the more distant 
low-order mean motion resonances may remain 
trapped during the growth of the protoplanet. For 
example, the resonant capture condition is most 
easily satisfied for the "distant" 3:2 and 2:1 mean 
motion resonances. As the protoplanet grows, the 
resonance strengthen increases and the libration 
time scale is reduced. Both the resonant prob- 
ability and the characteristic eccentricity of the 
resonant planetesimals increase. In the standard 
model with Tgrow = 10^ yr, the mass doubling time 
scale prior to the termination of growth is ~ 10'^ 
yr, which is comparable to the libration time scale 
in the mean motion resonances. Planetesimals 
that are loosely bound to the mean motion reso- 
nances have longer libration time scale than Tres- 
They are shaken by the rapid evolution of the 
protoplanet 's gravitational potential and cannot 
respond through adiabatic adjustments. This im- 
pulse leads to a late episode of planetesimal bom- 
bardment, which can introduce metallicity and 
structure diversity to the gas giant planets. The 
impact of this impulsive shake up is most pro- 
nounced in the rapid gas accretion models (see 
Figs. 8a and 8d), where the resonant capture 
becomes ineffective. 

The condition for resonant trapping is also more 
easily satisfied for the intermediate mass plan- 
etesimals, because their eccentricity damping and 
orbital migration time scales are relatively long. 
These tendencies are clearly evident in Figs. 3 
and 4. During the epoch of gas depletion, the 
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damping time scales Taero and Ttidc lengthen, which 
again provide favorable conditions for the capture 
of residual planetesimals into the mean motion res- 
onances. 

Finally, the excess density in the mean motion 
resonances is determined by the migration rate 
outside the resonances. Through orbital decay, 
planetesimals from the regions outside the first- 
born protoplanet congregate near its mean motion 
resonance. The enhancement of the local surface 
density promotes the formation of secondary gas 
giants. 

4. Summary and discussion 

In this paper, we investigate numerically the 
process and efficiency of planetesimal accretion 
onto a growing protoplanet. We consider the stage 
after the formation of the protoplanet core and 
during the accretion of its envelope. We use a 
physical Bondi formula and an ad hoc linear pre- 
scription to evaluate the gas accretion rate. The 
seed protoplanet is placed at 5 AU with initial 
mass 5.67 M^, accreting gas in a time scale Tgrow. 
The results of our numerical simulation and anal- 
ysis have implications for three issues: 
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Fig. 11. — Probability of survivers, colliders, 
ejectors(a > 100 AU) and Sun cashers (a < 1 AU) 
as functions of the planetesimals' mass m in the range 
lO^'' - 10^'' g at f = 2 X lO'^yr. The gas-accretion 
model is Bondi accretion. The model parameters are 
Tgrow = 10''' yr, protoplanet density p = 0.125 g cm~^. 



4.1. Suppression of planetesimal accretion 
during the onset of gas accretion. 

We first examine the initial growth of the pro- 
toplanet through gas accretion (transition from 
stages of embryo-growth (SI) to quasi-hydrostatic 
sedimentation (S2)). The accretion rate of the 
gas is suppressed by the bombardment of plan- 
etesimals, which generates heat needed to be re- 
distributed efficiently. In the vicinity of the core, 
there are two regions that supply the bullet plan- 
etesimals. Through its secular perturbation, the 
protoplanet induces planetesimals within an in- 
ner feeding zone to attain radial excursion, which 
crosses its orbit during each azimuthal conjunc- 
tion. Planetesimals within this band (~ l.S/iOp) 
engage in repeated close encounters. At 5AU, the 
local Keplerian speed is comparable to the surface 
escape speed of the protoplanet, so only a frac- 
tion of the close encounters will lead to physical 
collisions and the buildup of the core. 

The protoplanet also has an outer feeding zone 
at 1.5 — 3.5/iap. Due to the tidal perturbation of 
the host star and the protoplanet, planetesimals in 




Fig. 12. — Configuration of tidal drag perturbation 
on a planetesimal (see Appendix). 
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this region occasionally cross the orbit of the pro- 
toplanet. However, the frequency of such encoun- 
ters decreases with 5a and vanishes at the bound- 
ary of the feeding zone. In a gas frcc-cnvironmcnt, 
many synodic periods are needed for the planetes- 
imals in the outer feeding zone to be captured by 
the protoplanct. The presence of gas can lead to 
eccentricity damping, orbital migration, and dif- 
fusion of planetesimals from the outer to the in- 
ner feeding zone. The migration time scale is de- 
termined by the gravitational perturbation of the 
protoplanet as well as the efficiency of gas damp- 
ing. At the onset of the gas-accretion stage, the 
mass of the protoplanetary core is relatively low, 
so the excited eccentricities of the planetesimals 
arc relatively small. Consequently the migration 
rate from the outer to the inner feeding zone is 
slow. 

The low replenishment rate as well as the mod- 
est collision probability imply that the planetes- 
imal bombardment rate onto the core is much 
less frequent than that inferred from the efficient 
consumption of all planetesimals engulfed by the 
expanding feeding zone (as assumed in the early 
models of Pollack et al. 1996). The suppression of 
planetesimal collisions also lead to a decline in the 
energy dissipation rate and a cutoff in the replen- 
ishment of grains in the gaseous envelope of the 
protoplanet. The elimination of these bottlenecks 
for gas accretion shortens the growth time scale 
for proto-gas-giant planets. 

4.2. Chemical and structural diversity 
through late-stage bombardment 

A substantial fraction of the super-solar chem- 
ical composition in Jupiter resides in its enve- 
lope. We explore the possibility that the chem- 
ical and corc-cnvclopc structural diversity may be 
due to late-stage bombardment by residual plan- 
etesimals. Due to the secular perturbations of the 
protoplanct and gas drag, exterior planetesimals 
(with a > Op) in the outer feeding zone diffuse into 
the inner feeding zone and engage in close encoun- 
ters with the protoplanct. The interior planetes- 
imals (with a < Op) undergo orbital decay out of 
the feeding zone. Both effects lead to the clearing 
of the feeding zone and the formation of a plan- 
etesimal gap. 

In tenuous regions of the disk, protoplanet 
grows gradually and the orbits of both planetesi- 



mals and the massive embryos decay slowly. The 
migration of the external planetesimals and em- 
bryos may be stalled near the outer mean mo- 
tion resonances of the protoplanet. But during 
the growth of the protoplanet, its feeding zone 
expands and the mean motion resonances (with 
high p's) overlap. The planetesimals accumulated 
in the resonances become dynamically unstable. 
The colliding embryos may penetrate deeply into 
the protoplanet envelope and become part of its 
core. 

However, in the relatively dense inner regions of 
the disk, the gas-accretion rate of protoplanet and 
the orbital decay rates of the planetesimals are 
both relatively high. Embryos pass through the 
mean motion resonances without any significant 
perturbation. Only the intermediate-mass plan- 
etesimals can be captured onto the mean motion 
resonances of the protoplanet. When the proto- 
planet becomes sufficiently massive to destabilize 
the resonance-trapped planetesimals, the plan- 
etesimals will collide with the protoplanet. How- 
ever, as they do not have adequate mass to survive 
the passage through its envelope, they mostly sup- 
ply the metallicity in the protoplanet envelope. 

These two possible outcomes may account for 
the structural diversity between Jupiter and Sat- 
urn. Jupiter formed early in a relatively dense 
region just beyond the snow line. Intermediate- 
mass planetesimal bombardments occurred dur- 
ing the advanced stage of its progenitor's growth 
so that it has a relatively low core mass but a 
metal-rich envelope. In contrast, Saturn formed 
during the depletion epoch in the relatively tenu- 
ous outer regions of the disk. In this case, mas- 
sive embryos may be trapped in the mean motion 
resonances of the protoplanet when its mass is a 
few Mq. Their orbits become unstable when Sat- 
urn acquired most of its present-day mass. The 
late-stage bombardment by these massive embryos 
may have contributed to substantial core of Sat- 
urn. 

Generally the longer the time scale of gas ac- 
cretion, the more efficient a protoplanct accretes 
planetesimals. But when the time scale is compa- 
rable to the age of the gaseous disk (millions of 
years) , a runaway type of gas accretion model like 
Bondi accretion is preferred for the protoplanet 
to acquire more efficient and late-stage planetes- 
imal accretions (Figs. 8-9). In a solid disk with 
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surface density twice the minimum solar nebula, 
the solid mass colliding with a compact proto- 
planet is 6 '--^ 7 , and it reaches maximum when 
T £ [10^, 10^] yr in the Bondi gas-accretion model. 
This mass is comparable to the initial core mass 
and should be regarded as a lower limit. In our 
calculations, the radius is calculated according to 
equation (fTTj) . After the protoplanet acquired an 
atmosphere, the collision rate could be enhanced 
(Inaba & Ikoma 2003) and a higher solid accretion 
rate may be expected. 

The accretion rates of planetesimals with differ- 
ent masses are determined. The mass of individ- 
ual planetesimal that could be effectively accreted 
to the protoplanet lies in the range 10^° — 10^^ 
g, which corresponds to an embryo with radius of 
50 — 5000 km (Fig. 10). The relatively low mass 
planetesimals disintegrate in the protoplanet en- 
velope, whereas massive protoplanetary embryos 
may survive their passage through the envelope 
and become a part of the protoplanet core. 

4.3. The enhanced formation of multiple 
gas giant planetary systems 

We also show that due to the gas drag and pro- 
toplanet secular perturbations, the density of the 
planetesimal disk could be slightly increased out- 
side the orbit of the protoplanet, which will en- 
hance the subsequent formation of external proto- 
planet cores (Figs. 3-5). 

During their orbital decay, all planetesimals, 
with the exception of the most massive em- 
bryos, migrate sufhciently slowly that they be- 
come locked onto the mean motion resonances of 
the protoplanet. When the growth of the proto- 
planet is stalled, the resonant planetesimals re- 
siding outside the asymptotic feeding zone remain 
well separated from other mean motion resonances 
and are stably trapped in these resonances. The 
enhancement of the local surface density reduces 
the growth time scale of the embryos. The forma- 
tion of second generation proto-gas-giant planets 
is promoted. 

The emergence of multiple gas giants close to 
each other's mean motion resonances may also 
lead to dynamical instabilities (Zhou et al. 2007). 
The resulting dynamical interaction may lead to 
mergers, ejections, and breakup of the system. We 
suggest that this may be a promising avenue for 



the generation of the large eccentricity distribu- 
tion among the extra solar planets. 
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A. Perturbations under tidal drag 

Suppose the perturbing acceleration of tidal drag to a planetesimal has the form of equation (fT5 



/tidal = - ^' ^" =A(Vfc-Vg), (Al) 
Ttidal 

where A = l/rtidai- The velocity of the planetesimal and gas can be expressed in terms of the radial, 
azimuthal and normal components with unit vectors f,Tp and h, respectively (Fig. 12): 

Vfc = t;o(cosa f + sin a ■(/>), ^^^^ 
Vg — x(cos e V' + sin e h), 

where a is the angle from the radial to the velocity direction in the planetesimal orbit plane, and x = 
^^GMqI R{\ — 2rj{R)y^^. Suppose the gas is in a circular orbit, rj = 0, we obtain, 

^^„a(i±££^)i/2cos-V2 5. (A3) 
1 — e'' 

The definitions of the angles are also shown in Fig. 12. From spherical geometry (Adachi et al. 1976), 

sine = cos i/ cos (5, 

cose = sinicos(/ + ct;)/cosJ, (A4) 
cos (5 = [1 -sin2isinV + ^)]^^^ 

where / is the true anomaly of the planetesimal orbit. Thus the perturbing acceleration of tidal drag can be 
expressed as: 

/tidal = ^;^{esin/f + (1 + esin/)i/2[(i + esin/)i/2 _ cos z cos'^/^ 

+ (1 + esin/)i/2 sini cos(/ + uj) cos-3/2 Sh} (A5) 
EE Rf + T^j + Nh 

The evolution equations of the osculating elements under tidal drag obey (Murray & Dermott 1999): 

t =;^7^[i?esin/ + r(l + esin/)] 

= + e2 ^2ecos/- (1 + e sin /)3/2 cosi[l - sin^ i sin2(/ + tj)]-3/4} 

rfe ^ [j^ sin / + r(cos / + )] 

at na L j ' \ j ' l+ecos/''-l (A6) 

= 2A(e + cos/) + A(l + ecos/)-3/2(2cos/ + ecos2/ + e)cosi[l -sin2isin(/ + w)]-3/4 

di _ yi-e^ cos(/+a)) 
dt na 1+e cos / 

= Asini(l + ecos/)-i/2[l - sin^ i sin(/ + uj)]-^/'^ cos2(/ + Q). 



We average over time t in one period of Keperian motion according to : 

TJo 2ttJo (1-ecos/) 

This gives 



<F>=- I Fdt^ — I j-}——l^df. (A7) 



i<rf|> =A[{l-§e^-^t^ + ^t^sm^uj) + o{e^,i^)], (A8) 
T < f > = 41(1 - i^' + ie^sin^o; + ^^') + oie',^')]. 

We further average over one period of du/dt to eliminate the dependence of lo, and recall A = — 1/rtidai, 

which finally yields 

i<t> =-^l(l-i^'-i^') + «(e^^^)], (A9) 
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